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The 20 1 1 Tohoku earthquake (M w = 9.1) highlighted previously unobserved features for megathrust events, 
such as the large slip in a relatively limited area and the shallow rupture propagation. We use a Finite 
Element Model (FEM), taking into account the 3D geometrical and structural complexities up to the trench 
zone, and perform a joint inversion of tsunami and geodetic data to retrieve the earthquake slip distribution. 
We obtain a close spatial correlation between the main deep slip patch and the local seismic velocity 
anomalies, and large shallow slip extending also to the North coherently with a seismically observed 
low-frequency radiation. These observations suggest that the friction controlled the rupture, initially 
confining the deeper rupture and then driving its propagation up to the trench, where it spreads laterally. 
These findings are relevant to earthquake and tsunami hazard assessment because they may help to detect 
regions likely prone to rupture along the megathrust, and to constrain the probability of high slip near the 
trench. Our estimate of —40 m slip value around the JFAST (Japan Trench Fast Drilling Project) drilling 
zone contributes to constrain the dynamic shear stress and friction coefficient of the fault obtained by 
temperature measurements to —0.68 MPa and —0.10, respectively. 

On March ll lh 2011 one of the largest earthquakes ever recorded occurred at the subduction interface 
between the Pacific and the Okhotsk plates and struck the Tohoku region in Japan (Fig. 1). This M w 9.1 
earthquake, located at 142.68°E 38.19°N\ generated a tsunami that devastated the Japanese coasts, 
including towns and important infrastructures such as the Sendai airport and the Fukushima nuclear power 
plant causing more than 16,000 fatalities 2 . The Tohoku earthquake is also the best observed ever megathrust event 
and consequently it has been investigated by modelling the unprecedented high-quality data set recorded by the 
Japanese dense seismological, geodetic and marine observational networks. The numerous studies published in 
recent literature are based on different kind of data and methodologies, including teleseismic 3 , strong motion 4 , 
geodetic", tsunami waveforms 2,8 , and joint inversions 911 which were performed to investigate the earthquake 
rupture process. 

The resulting source models share two common features of the coseismic rupture, stimulating further investi- 
gations to explore the physical processes controlling the genesis and the impact of megathrust events. The first 
feature is that the overall Tohoku rupture area is mainly concentrated in a relatively small portion of the plate 
interface and the retrieved peak slip values range between 30 and 60 m. This long-wavelength feature is common 
to most of the slip models obtained using different kind of data. At the same time, differences in terms of 
maximum slip value or number of slip patches can be observed among models, due to the data resolution and 
fault parameterization 4,7,10,12,13 . A rather small rupture area characterized by very large slip is quite unusual for 
great earthquakes such as the 2004, M w = 9.2, Sumatra 14 ' 15 and the 2010, M w = 8.8, Maule 16 megathrust events. 
Some authors observed also a first order correlation between the coseismic slip patch and the positive seismic 
velocity anomaly at the subduction interface 1720 . The relatively concentrated Tohoku rupture area maybe related 
to the lithosphere structure and the consequent heterogeneous pattern of pre-stress as well as to the fault frictional 
properties of the plate interface that could promote or inhibit the rupture propagation. The relative variations in 
shear-wave and bulk-sound speed detected in the coseismic slip area may reflect mechanical heterogeneities of the 
subduction interface, which may have acted as asperities allowing this event to build up large slip in the near- 
trench zone 21 . 
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Figure 1 | Location Map. (a) Red star indicates the epicentre position. Red and white "beach ball" represents the focal mechanism of 2011 Tohoku 
earthquake. Yellow triangles indicate the DART stations used in the inversion; (b) Cyan circles indicate GPS stations onshore, magenta circles indicate the 
geodetic seafloor observation sites, yellow triangles indicate the bottom pressure sensors and GPS-buoys (Table S4 in Supplementary Information) . White 
arrow indicates the approximate convergence direction of the Pacific plate (estimated velocity of 9.2 cm/yr). Maps are created using GMT (Generic 
Mapping Tools, http://gmt.soest.hawaii.edu/) software. 



The second intriguing feature of the 2011 Tohoku earthquake is 
the shallow rupture propagation up to the trench 2,4,10,13 . The large 
(>15 m) and shallow coseismic slip is at odds with the quite com- 
mon view on the coseismic behaviour of the shallow portion of 
megathrusts, where aseismic slip, lack of strain accumulation and 
low coupling are expected to be dominant 22 . Recent results from the 
Japan Trench Fast Drilling Project (JFAST 23 " 25 ) pointed out the pres- 
ence of smectite-rich weak clay, up-dip from the hypocentre in the 
very shallow portion of the subduction plate boundary. Thus, this 
distinctive feature of the Tohoku earthquake might have been con- 
trolled by very low friction 26 on a relatively thin fault zone rich of clay 
sediments 24,25,27 . 

Both these features, contributing to confine a vertical seafloor 
displacement larger than 10 m in a relatively small area 2,10 , may have 
increased the tsunamigenic potential of the 2011 earthquake, gen- 
erating tsunami waves higher than 10 m and runups larger than 
30 m along the Iwate prefecture coasts 28 . 

Accounting for detailed information on the geometry and struc- 
ture of the subducting plate is necessary in order to reduce the epi- 
stemic uncertainties related to the modelling, to get a robust image of 
the coseismic slip distribution, to further constrain its spatial exten- 
sion with respect to previous studies, and to focus on the shallow 
near-trench portion of the megathrust. We achieve this by carrying 
out, for the first time to our knowledge, a joint inversion of tsunami 
and both inland and seafloor geodetic data, constraining the slip 
distribution by means of a 3D Finite Element Model (FEM) of the 
subduction zone. We characterize each single element of the FEM 
grid by the 3D elastic structure inferred from seismic tomography 29 , 
whereas the shallow near-trench portion of the megathrust, rich of 
sediments, is modelled by assuming a much more compliant material 
with respect to the surrounding medium (see section Methods). Our 
modelling approach enables to account for the effects of the elastic 



contrasts both within the crust and between the crust and the mantle 
and, furthermore, for mimicking the contrast between the crust and 
the sediments at the trench (i.e. clay-rich sediments), which may play 
a key role in controlling the coseismic slip distribution. 

Aim of this work is to image the slip distribution adopting a more 
realistic representation of the Green's functions through the 3D FEM 
model. In particular, we discuss whether and how the rupture extent 
and its propagation near the trench have been controlled by the 
regional and fault-zone structural heterogeneities. We also relate 
the shallow slip propagation to the recent JFAST results, drawing 
some possible relationships with the observed seismic radiation. 

Results 

The 3D geometry of the FEM (cf. section Methods for details) is built 
taking into account the main features of the subduction zone includ- 
ing the Tohoku region, such as the slab geometry and the topo- 
graphy/bathymetry. A FEM mesh of about 280,000 8-nodes brick 
elements is then created (Fig. 2). The 3D elastic structure of each 
element is constrained by the 3D P-wave and S-wave velocities (Vp 
and Vs, respectively) of a regional tomographic model 29 ; the resulting 
crust, slab and mantle are highly heterogeneous, with rigidity values 
ranging from —20 GPa in the crust, to —50-60 GPa in the slab. The 
elastic structure of the shallow part of the megathrust (i.e. the prism 
above the subduction fault) is characterized by independent data as 
imaged by seismic reflection surveys 30,31 . The prism volume is limited 
by the trench to the East and by a steep listric plane to the West 30 and 
is considered as made of unconsolidated clay-rich sediments and 
volcanic materials, for which we assume uniform elastic constants 
(rigidity = 5GPa and Poisson's ratio = 0.33). The chosen elastic 
values are compatible with those attributable to weak-clay (smectite) 
observed by the JFAST drilling project 23 " 25 . 
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Figure 2 | FEM model, (a) Total domain of the 3D FEM model (—2900 X 2500 km 2 ). Yellow dots represent the free surface grid nodes, (b) Central part of 
the 3D FEM model, including the active fault, viewed from SW. The light blue line is the trench and the red line is the section of the active fault, (c) The 
inset shows a zoom of the interface between trench and the uppermost part of the fault. The shaded element edges lines on the top identify the 
prism extent. Details about the elastic layering of the model can be found in the Supplementary Information, Figure SI. Panel (a) is created using GMT 
software, panels (b,c) using AMIRA (http://www.vsg3d.com/amira/overview, Date of access: 1 1/06/201 1) software. 



The FEM model is then used to compute a set of Green's Functions 
which are combined with a robust joint inversion scheme adopted in 
several previous papers to retrieve the slip distribution 10,14,15,3233 . The 
inverse problem is solved by using the method of Green's functions 
superposition and the Heat Bath algorithm, a particular implementa- 
tion of the Simulated Annealing technique 34 . The slip distribution is 
obtained through a global search technique 35 . Since the solution of 
this inverse problem is intrinsically non-unique 36 , and to account for 
possible modelling uncertainties, we compute and show the average 
slip model that we consider more representative than the best fitting 
solution 16 . The average slip model is calculated as a weighted mean of 
a subset of the explored models. We assess the dispersion of the 
model parameters around their average values by performing an a- 
posteriori analysis of the explored models ensemble so that the error 
associated to each parameter is its weighted standard devi- 
ation 10,14,15,32,33 . In this analysis the coefficient of variation (i.e. the 
standard deviation divided by the average slip value) is also 



computed for each subfault parameter (further details in Methods 
section). 

The 2011 Tohoku earthquake slip distribution. The slip 
distribution of the 2011 Tohoku earthquake (Figure 3a) extends 
mainly along strike from ~36.5°N to ~39.5°N and up-dip from 
the hypocentre in the SE direction (from ~142°E to ~144°E). The 
portion of megathrust with slip values greater than —10 m (i.e. 
—20% of the maximum peak of slip) corresponds to an area of 
about 350 X 200 km 2 . Amplitudes ranging between 25 m and 
30 m characterize the slip distribution in the earthquake 
nucleation zone. Slip amplitudes increase moving away from the 
hypocentre towards the trench; slip direction is consistent with 
both the relative convergence of the Pacific and the Okhotsk plates 
and the centroid moment tensor solution (GCMT). Rake angle 
averaged over the whole slipping area is —88°. The rupture reaches 
the shallow and less locked part of the megathrust with slip values 
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Figure 3 | Slip distribution. Slip distribution for the 2011 Tohoku-oki earthquake obtained from the joint inversion of tsunami and geodetic data, (a) 
Orange arrows represent the slip direction (rake, Table S3 in Supplementary Information) . Thin dashed black contours above the fault plane indicate the 
interseismic coupling 48 (from 10% to 100%, at 10% intervals) along the megathrust. Black arrow indicates the approximate convergence direction of the 
Pacific plate (estimated velocity of 9.2 cm/yr). Red star as of Figure 1. (b) Green contour lines (intervals are at 0.25, 0.50, 1.00, and 1.50 m) and green star 
indicate the slip distribution and the epicenter position of the foreshock 45 occurred on 9 March 201 1, respectively; magenta dashed rectangle represents 
approximately the rupture area of the 1896 Meiji-Sanriku earthquake 2 ; yellow coloured region approximately indicates the zone of coseismic high 
frequency radiation; black dot indicates the JFAST ocean drilling site C0019 23 ; (c) coefficients of variation associated to the average slip model of the 2011 
Tohoku-oki earthquake resulting from the joint inversion. Black contour lines (lOmeters interval) indicate the slip distribution of Tohoku earthquake. 
Maps are created using GMT software. 
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greater than 30 m, up to ~45 m around ~143.5°E, 38°N. The 
inferred slip distribution shown in Figure 3b displays a narrow 
patch of slip located to the NE (from ~38.5°N to ~39.5°N) in the 
shallow part of the subduction interface roughly corresponding to 
the area where the 1896 Meiji-Sanriku earthquake occurred 37 . In 
addition, at around 35 km depth a relatively low slip feature 
extends from the hypocentre zone to the latitudes of Fukushima 
and Ibaraki prefectures. This feature, generally imaged by using 
strong-motion and seismic data 413 , approximately corresponds to a 
source of high frequency radiation detected by means of back- 
projection analysis 12 and amplitude source location method 38 . Even 
though our inverted distribution represents a long- wavelength image 
of the coseismic slip pattern, the complex model adopted and the 
inversion technique allow sensing this peculiarity of the Tohoku 
earthquake even using only tsunami and geodetic data. We 
recognize, however, that the pattern of the variation coefficient 
(Figure 3c) indicates that the slip distribution is particularly well 
resolved (variation coefficient lower than —0.4) in the portion of 
megathrust with slip values greater than 10 m. 

The estimated total seismic moment, computed by taking into 
account the elastic parameters used into the FEM model (see 
Methods section) is M 0 = 5.72 X 10 22 Nm, corresponding to a mag- 
nitude M w = 9.1. 

The inferred slip model (Fig. 3a) yields a general good fit to tsu- 
nami (RMS = 0.50 m) and geodetic (RMS = 0.38 m, 0.19 m, and 
0.14 m for East- West, North-South, and vertical components, 
respectively) data (Fig. 4). It is noteworthy that the tsunami signals 
measured by the GPS-buoys off Iwate (G802, G804, G807) and by the 
bottom pressure sensors located very close to the source are extre- 
mely well reproduced, in particular their short-period wavelengths. 
This feature is probably tied up to the narrow patch of slip located in 
the shallow part of the megathrust from ~38.5°N to ~39.5°N. 
However, a slight phase shifting is observed at few tsunami stations 
(e.g. G801) and it may be due to our assumption of a simplified 
circular rupture front (propagating from the hypocenter with a velo- 
city of 1 km/s 39,40 ) while the actual rupture history of Tohoku 
event 4,41 is likely more complex. 

We also notice, as already observed in previous studies 2,11 , that the 
tsunami amplitude is amplified by the contribution of the horizontal 
seafloor deformation near the trench axis (due to the presence of 
steep bathymetric slopes 42 , further details in section Methods). This 
effect is particularly evident for the stations located near or above the 
source region (TM1, TM2, P02, P06, G802, G804, G807) where the 
wave amplitudes are smaller by —20-40% if the horizontal deforma- 
tion contribution is not included (Supplementary Figure S3). 

The present slip model confirms at the first order the rupture 
pattern retrieved in our previous study 10 as well as in other 
works 2,4,5,7,9,1113,43,44 , i.e. the large slip in a relatively concentrated 
rupture area and the shallow rupture propagation. In particular, a 
shallow patch of slip in the northern portion of the megathrust is 
clearly distinguishable, a feature also observed using seismic 
data 4,12,13,44 and that will be discussed in a later section. In addition, 
the present slip model highlights a new interesting feature, thanks to 
the employment here of an enlarged geodetic and tsunami data set 
(particularly above the source, e.g. P02 and P06 sensors), and of a 
FEM model honouring the structural heterogeneities. Indeed, we 
observe a relative minimum of slip (<10 m) close to 38.5°N, 
143°E (Figure 3b) that overlaps with the main rupture area of the 
foreshock (M w = 7.4) occurred on 9 March 201 1 45 . This distinctive 
feature may be attributed to a low pre-stress level due to the fore- 
shock occurrence; actually, a negative stress change has been esti- 
mated in the portion of megathrust slipped during the foreshock 45 . 

Rupture extent, seismic velocity anomalies and frictional pro- 
perties. The mechanism of generation, propagation and arrest of 
an earthquake may depend on different factors such as 3D 



structural heterogeneities of the medium, the frictional properties 
of the source zone, and the stress load. The seismic velocity 
anomalies at the fault-scale could represent a reasonable proxy to 
describe the variability of the above-mentioned features 18,21,45 . 
Moreover, some authors 18 observed that in the rupture area of the 
Tohoku earthquake the seismic velocity anomalies are correlated 
with the slip distribution obtained by inverting only inland GPS 
data 47 . Here, we compare our slip distribution with the high- 
resolution tomographic model at the subduction interface 29 in 
order to better understand how much the Tohoku rupture has 
been controlled by structural heterogeneities (i.e. seismic velocity 
anomalies) at the fault-scale. 

Three distinct zones, remarkably evident for Vp anomalies distri- 
bution, are present along the seismogenic portion of the slab inter- 
face, in a well-resolved region of the tomographic model 29 . We 
observe in Figure 5a a Positive Velocity Anomaly zone (hereinafter 
PVA) centred at ~38°N, ~142.5°E, whereas two Negative Velocity 
Anomaly zones (hereinafter NVA) are located just to the North and 
to the South of the PVA (centred at ~39.5°N, ~143°E and ~37°N, 
~ 142. 5°E, respectively). A very strong spatial correlation is observed 
between the slip distribution (slip > 10 m) and the PVA present in 
the Vp model; the rupture area almost completely overlaps the PVA 
zone (>70% of the slipping area) and borders the NVA zones. Such a 
correlation holds for the Vs anomalies as well (Fig. 5b). This close 
correspondence suggests that the Tohoku rupture might have been 
efficiently controlled by the observed velocity anomalies, perhaps 
associated to the variation of the frictional properties of the 
Tohoku megathrust zone. 

Comparison with the map of interseismic locking 48 at seismogenic 
depth (>~20 km, where the coupling estimation is well resolved) 
shows that the coupling degree is very high (>~80%) in the PVA 
zone (Figs. 5a,c), while it decreases (<~50%) moving toward both 
the NVA zones. In principle, zones with relatively low coupling are 
assumed to be creeping and having a velocity strengthening beha- 
viour, while locked zones, where the earthquakes generally occur, are 
assumed to have a velocity weakening behaviour 49-51 . Hence, we infer 
that the NVA zones, characterized by relatively low coupling, and the 
PVA zone, that is the most locked, exhibit a velocity strengthening 
and velocity weakening behaviour, respectively. This interpretation 
is supported by the observation that past large interplate earthquakes 
(M > 7 occurred since 1900 along the Tohoku region) are distributed 
preferentially inside the areas characterized by PVA 18 as shown in 
Figures 5a,b. 

In summary, the absence of large past earthquakes, the low coup- 
ling, and the relatively confined spatial extent of the Tohoku earth- 
quake rupture area (mainly restricted to the PVA region) considered 
together suggest a possibility that the NVA may have acted as velocity 
strengthening zones. This hypothesis accounts for the unexpected 
and very high slip in a relatively small area, conflicting with expected 
values from empirical earthquake scaling laws 52 . The fact that a rela- 
tively small percentage of the rupture area (<30%) lies within the 
NVA zones is not at odds with this interpretation. Such a feature is 
not unusual for great earthquakes such as the one occurred in 
Tohoku. Indeed, a recent study has demonstrated that under specific 
conditions great earthquakes ruptures might propagate, even par- 
tially as in the present case, through velocity strengthening zones 51 . 

A further analysis of the seismic velocity anomalies provides 
insights about the possible fluids content along the subduction inter- 
face nearby the Tohoku hypocenter. Actually, the earthquake nuc- 
leation may be also related to the distribution of fluids around the 
plate boundary 53 . The earthquake nucleated in a strongly coupled 
zone that is characterized by both PVA and Positive Vp/Vs Ratio 
Anomaly (PVRA from here on, Fig. 5d). The PVRA may indicate 
presence of fluids in brittle rocks 54 . As a consequence, excess pore 
pressure might have decreased the strength on the fault plane and 
caused the main shock occurrence. 
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Figure 4 | Comparison between observed and predicted data sets. Comparison between observed (black) and predicted (a) horizontal and (b) vertical 
displacements at GPS (red) and geodetic seafloor observation sites (magenta), (c) Comparison between the observed (black) and predicted (red) tsunami 
waveforms. Panels (a,b) are created using GMT software; panel (c) using MATLAB (http://www.mathworks.it/products/matlab/, Date of access: 11/06/2011). 
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Figure 5 | Comparison between Slip and Velocity anomaly distributions. 

(a) Vp and (b) Vs anomaly distributions 2 * on the subducting plate 



interface. Black contour lines (10 m interval) indicate the slip distribution 
of Tohoku earthquake. Open green circles indicate the large (M > 7) 
earthquakes occurred in the Tohoku earthquake region since 1900. Orange 
dashed contour lines ( 1 8 km interval) indicate the depth of plate interface. 
Red star as of Figure 1. (c) Slip distribution for the 2011 Tohoku-oki 
earthquake. Red star and thin dashed black contours above the fault plane 
as of Figure 3a; (d) Vp/Vs ratio anomaly distribution on the subducting 
plate interface. Black contour lines (10 m interval) indicate the slip 
distribution of Tohoku earthquake. Red star as of Figure 1 . Maps are 
created using GMT software. 

In the light of all of the above, we suggest that the velocity strength- 
ening zones may have initially acted as barriers to rupture propaga- 
tion, laterally confining it to the relatively small area around the 
earthquake nucleation point. Subsequently, due to the coseismic 
weakening and the local very heterogeneous frictional condi- 
tions 24,25,27 , the rupture has been driven towards the trench where it 
further propagated and spread laterally along strike in the low fric- 
tion shallow part of the subduction zone, as shown by the slip model 
(Fig. 3a). The shallow slip propagation will be more thoroughly 
analysed in the next section. 

Shallow slip propagation. The slip distribution retrieved in the 
present study and most of the published source models 2 ' 4,9 ' 10 ' 13,55 
highlight large slip values (>15 m) in the shallower part of the 
subducting interface (depth <~ 15 km) from ~37°N to ~39.5°N, 
with peak values greater than —35 m nearby 143. 5°E, 38°N. Since the 
shallowest part of megathrusts is commonly thought to be 
characterized by a rate strengthening behaviour 22 and dominated 
by aseismic slip, lack of strain accumulation and low coupling, the 
large shallow slip occurred during the Tohoku earthquake came as a 
surprise. 

Numerical modelling of the dynamic rupture 56 and laboratory 
experiments on clays 27 hypothesized dynamic weakening of the fault 
at high slip rate as the possible reason for the large slip near the 
trench. This hypothesis is corroborated by the results of JFAST pro- 
ject 23 25 that conducted a drilling survey near the Japan trench (site 
COO 19, Fig. 3b). In particular, dynamic shear stress and friction 
coefficient were determined by the temperature measurements at 
the plate boundary 24 and by laboratory friction experiments 25 under 
conditions of both low permeability and high slip velocity. The 
JFAST results as a whole indicate low shear stress, low friction on 
the fault, and low permeability conditions of smectite-rich clays 
within the decollement, thus supporting the coseismic weakening 
of the fault. The estimation of the friction coefficient from temper- 
ature has been obtained by assigning a value for slip, which is one of 
the greatest uncertainties in the estimates of shear stress and friction. 
A fault slip of 50 m is assumed 24 at the site C0019, in agreement with 
some previously published slip models. In this area our model indi- 
cates instead a lower slip value of ~40 m (Fig. 3a), which would 
imply average shear stress and friction coefficient of —0.68 MPa 
and —0.10 respectively (Patrick Fulton, personal communication). 
These values are slightly higher with respect to estimates using 50 m 
of slip but still in the range of values obtained experimentally for slip 
greater than — 35 m under different permeability conditions 25 . 

The slip distribution found in the present study features a rupture 
pattern further extending to the northern portion of the Tohoku 
region, in a narrow and shallow stripe from — 38.5°N to — 39.5°N, 
with slip values up to 26 m (Fig. 3b). The northern patch of slip 
roughly overlaps the estimated rupture area of the M w = 8.0-8.2 57 
1896 Meiji-Sanriku earthquake 2 ' 37 . The maximum tsunami runup 
heights caused by the Tohoku event along the Iwate prefecture coasts 
and those caused by the 1896 Meiji-Sanriku tsunami earthquake 58 are 
comparable 28 . Some authors 2,8 examined possible similarities 
between the causative source of this shallow past tsunamigenic event 
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and the shallow part of the Tohoku rupture. The 1896 event has been 
characterized by a very weak high-frequency seismic radiation and a 
long rupture duration 58 , and it has been classified as a tsunami earth- 
quake. A slow component in the source of the 201 1 Tohoku event is 
not detected by using long-period spheroidal modes of the Earth 59 , 
but this could be likely due to the lack of sufficient resolution at a 
local scale. Contrarily, different images of the radiated energy 
obtained in recent studies about the Tohoku earthquake are compat- 
ible with the hypothesis of slow rupture for the northern patch of slip. 
Indeed, in this zone a weak high-frequency seismic radiation feature 
has been observed 43 , whereas a strong low-frequency seismic radi- 
ation contribution has been identified to the North of ~38.5°N 12 . 
The simultaneous presence of strong low- frequency and weak high- 
frequency radiation in the area to the north of ~38.5°N may indicate 
a relatively slow rupture 60 . This hypothesis is also consistent with the 
low slip rate in the same region 44 . Therefore, a shallow slow rupture 
with associated very large slip, in principle, is compatible with a 
tsunami earthquake mechanism. Nevertheless, some authors 6 claim 
that a submarine landslide is required in place of this patch of slip. 
Their slip model, derived only by geodetic data inversion, which 
clearly may have only a local resolving power on the slip distribution, 
indeed fails in predicting the tsunami signals at GPS-buoys off Iwate 
prefecture (i.e. GPS-buoys 802, 804, and 807) and therefore they need 
to invoke an additional tsunami generation mechanism to reproduce 
these data. The impulsive and short-period tsunami signal measured 
by the above-mentioned GPS-buoys and both TM1 and TM2 sensors 
are probably due to the short wavelength offshore seafloor deforma- 
tion generated by the large and shallow slip in the Tohoku region 2 . 
Therefore, we observe that the tsunami signals observed at GPS- 
buoys off Iwate (Figure 4c) are in principle fully compatible with a 
narrow and shallow patch of slip as the one in Figure 3, derived by 
our joint inversion. 

In conclusion, the coseismic slip is overall coherent with the 
observed seismic radiation and it is sufficient to reproduce the 
short-period tsunami signals without invoking additional tsunami 
sources. Even if this reasoning does not rule out the hypothesis of a 
landslide, it shows at least that such a contribution to the tsunami 
source is not strictly necessary to explain the observed data, as the 
differences with synthetic data are likely to be well within epistemic 
uncertainties. In order to be confirmed, the hypothesis of such a huge 
mass movement would need additional supporting evidence, e.g. a 
dedicated bathymetric survey. 

Discussion 

We derived the slip distribution of Tohoku earthquake by perform- 
ing a joint inversion of tsunami and geodetic data combined with 3D 
FEM technique. This novel approach based on high-quality data, and 
the modelling of geometrical and structural heterogeneities of the 
source zone permitted us to obtain an accurate image of the slip, 
particularly near the Japan trench. Our slip distribution is character- 
ized by very large slip values (up to —45 m) in a relatively concen- 
trated rupture area that extends up to the trench with slip values 
greater than 30 m. The features of the retrieved slip model are com- 
patible with the sources of high- and low-frequency seismic radiation 
of the Tohoku earthquake identified by back-projection 12 and source 
location methods 38 . 

Moreover, we investigated whether the slip distribution might 
have been controlled by structural or frictional heterogeneities at 
the fault-scale. The slip image shows an unprecedented spatial cor- 
relation with the positive seismic velocity anomalies, with the rupture 
area overlapping by more than 70% the positive anomaly zone and 
bordering the negative anomalies. We suggest that the negative 
anomaly zones have acted as velocity strengthening barriers for the 
rupture propagation. This hypothesis accounts for the unusual shal- 
low and very large slip in a relatively small area, which consequently 
may have increased the tsunamigenic potential of the earthquake. 



Finally, we examined the shallow rupture propagation of the earth- 
quake. The large and shallow slip values we found in the central part 
of the fault are compatible with the results of JFAST drilling sur- 
vey 2325 , even if our slip value (—40 m) indicates slightly higher 
dynamic shear stress and friction coefficient. Consistently with the 
results of strong-motion and seismic data inversions, we also found 
a shallow slip patch along the fault portion roughly corresponding to 
the rupture area of the 1896 Meiji-Sanriku tsunami earthquake. This 
patch of slip is reconcilable with some images retrieved by back-pro- 
jection technique corroborating the hypothesis of a slow rupture com- 
ponent and, implicitly, a tsunami earthquake mechanism. Our results 
also show that a seismic source is sufficient to explain the Tohoku 
earthquake rupture process and to reproduce the ensuing short-period 
tsunami signals without invoking additional tsunami sources. 

In conclusion, our findings encourage further studies to verify the 
correlation between the velocity anomalies and the megathrust rup- 
tures, as well as the trench structure and material properties, in other 
subduction zones. This kind of studies can be of utmost importance 
for long-term earthquake hazard assessments, because they may help 
to detect regions likely prone to rupture along the megathrust, and 
even more for tsunami hazard assessments, if they can help con- 
straining the probability of high slip in the trench region. 

Methods 

Tsunami data. Tsunami waves generated by the Tohoku earthquake have been 
recorded by many instruments positioned near the Japanese coasts (GPS-buoys, 
coastal wave gauges, bottom pressure sensors), above the source area (bottom 
pressure sensors) and in the open ocean (DART buoys). 

The original signals include the tide and in some cases the surface waves generated 
by the earthquake. Tides are removed by using a procedure based on robust LOWESS. 
Further details about tsunami data are in Supplementary Information. 

Geodetic data. The stations of the GPS Earth Observation Network (GEONET) of 
Japan have recorded the coseismic deformations associated with the Tohoku 
earthquake. GPS data used in this study are distributed mainly along the Honshu 
Island and the southernmost part of the Hokkaido Island where a distinct coseismic 
deformation has been measured. The Geospatial Information Authority (GSI) of 
Japan has provided all original data and their processing is described in a previous 
work published by the same authors 10 . 

The large coseismic deformation has been estimated also on the seafloor by some 
geodetic observatories operated by Japan Coast Guard 61 and by Tohoku 
University 62,63 by a GPS/acoustic combination technique. The coseismic seafloor 
deformations have been measured as difference between the positions observed 
before and after the 2011 event exploiting the acoustic transponder technology. 
Furthermore, additional vertical seafloor geodetic measurements have been obtained 
by some ocean bottom pressure gauges positioned above the seismic source 7,62 . 

FEM. The free surface of the FEM is determined by means of a 250 m spatial 
resolution topographic/bathymetric digital model provided by the Hydrographic and 
Oceanographic Department (HOD) of the Japan Coast Guard. The data cover the 
whole top area of the FEM domain, about 2900 km X 2500 km (from 125°E to 163°E 
and from 15°N to 52 C N, Fig. 2). The geometry of the megathrust is constrained by the 
slab model for subduction zones (Slabl.0 64 ), available at http://earthquake.usgs.gov/ 
research/data/slab/#models (Date of access: 11/06/2011). The obtained interface 
extends for 550 km from 35.4 N to 40.4°N, latitudes of intersection with the trench, 
and it follows the curvature depicted by the slab down to about 70 km depth, which is 
the supposed fault's bottom (Fig. 2). Along the Japan trench the topographic/ 
bathymetric and the slab surfaces have a variable vertical distance ranging between 
1 km and 3.5 km. The discrepancy at the trench is due to the different origin of the 
two 3D digital surface models and to the degraded slab surface resolution at the 
trench. The slab, whose thickness is 45-70 km, penetrates down to 600 km depth 
with a dip angle of A5° 17a9 . 

The spatial resolution of the grid is 3 km close to the trench and 6 km in the 
remaining fault, while the element dimension increases toward the lateral and bottom 
edges of the computational domain. The fault plane is subdivided into patches of 
variable size: 24 km X 14 km (length X width) close to the trench (up to —15 km 
depth), 24 km X 24 km in the central part up to —40 km depth (corresponding to 
more than 2/3 of the fault surface) and 35 km X 35 km in the deeper part. Null 
displacements are assigned to the bottom boundary of the numerical domain. The 3D 
elastic structure is constrained by the 3D Vp and Vs local earthquake tomography 29 , 
where the whole Honshu arc (from the Japan trench to the backarc in the Japan Sea) is 
mapped by inverting a high-quality data set of P- and S-wave arrivals, both from 
earthquakes prior to the Tohoku event and occurring under the land area and seismic 
events beneath the Pacific Ocean and Japan Sea. The suboceanic events help to 
constrain the 3D velocity model under the Pacific Ocean more reliably. The absolute 
seismic velocities are transformed into elastic parameters using 
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and 

V S <» 

where {)., fi) are the Lame's constants (medium considered isotropic) and p is the 
density. The density is assigned to each layer defined 65 66 in the FEM model 
(Supplementary Information, Fig. Si, Table SI). Note that each element of the grid is 
characterized independently based on the 3D seismic velocities 67 , without layering 
approximation that is referred to the density only. 

The 10 km vertical resolution of the tomographic model is too coarse to allow 
constraining the seismic velocities at shallow depths close to the trench, where more 
compliant material with respect to the surroundings is likely present 29 . Therefore we 
use independent data from seismic reflection surveys 303 ' to properly describe the 
volume above the fault and across the trench, considering the shallow prism volume 
as made of sedimentary (such as clay-rich sediments) and volcanic materials, for 
which we assume uniform elastic constants 30 ' 31 . 

Geodetic Green's functions. Green's functions at GPS stations and seafloor 
observation points are generated by means of the FEM. The active fault is subdivided 
into 398 subfaults (Supplementary Information, Fig. S2, Table S2) composed of pairs 
of coincident nodes. Each Green's function is computed by imposing kinematic 
constraints at the nodes of the FEM grid pertinent to the considered subfault. Since 
each subfault is made of several dislocating nodes, the slip may reach very shallow 
depths, close to the trench (~1 km). The coincident nodes of the active subfault 
dislocate of a relative fixed amount Au along the rake direction, while moving 
together in the perpendicular directions. The remaining pairs of coincident fault 
nodes are constrained to move together of the same amount (i.e., no dislocation 
among coincident nodes). The displacement is not distributed as — Au/2 and +Au/2 
at both sides of the subfault, but the amount at each side depends on the subfault dip, 
the distance from the free surface, and the local elastic parameters contrast. FEM 
computations are carried out using the commercial software Abaqus, version 6.9 6S . 

Tsunami Green's functions. The initial condition for tsunami propagation is 
generally represented by the water displacement associated to the vertical seafloor 
deformation. However, in the regions of steep bathymetric slopes the contribution of 
the seafloor coseismic horizontal deformation may be not negligible with respect to 
the vertical one 11,42 . Thus, the tsunami initial condition is computed as 

dH dH 

u = u z + u x — — \-Uy— (3) 
ox oy 

where u Xi u y , and u z are the components of the coseismic displacement in the east, 
north and vertical direction respectively, H is the bathymetry depth (positive 
downward) and u is the final tsunami initial condition. The components u x , u yi and u z 
are numerically computed by Abaqus as well. The tsunami Green's functions for each 
subfault then are modelled by means of NEO WA VE 69,70 , a nonlinear dispersive wave 
model for tsunami propagation. The code numerically solves the nonlinear shallow 
water equations by using a semi-implicit finite difference technique on a staggered 
grid. In order to take into account the weakly dispersive behaviour of tsunami 
propagation, in NEOWAVE the equations include a non -hydro static pressure term 
and a vertical momentum equation. For tsunami modelling we use a computational 
domain with a spatial resolution of larc-min, and the bathymetric digital model 
provided by HOD opportunely resampled. 

Inversion. We follow a joint inversion scheme whose reliability has been already 
implemented in several previous papers 10,14,15 ' 32,33 - 71 to retrieve the coseismic slip 
distribution of the 2011 Tohoku-oki earthquake (further details in Supplementary 
Information). Two different cost functions for tsunami and geodetic data are used. 
For tsunami data we use a function that results to be sensitive in matching both phases 
and amplitudes of the time series 72 . Instead, a standard L2 norm is used for the 
geodetic data set. 

Previous studies on the Tohoku earthquake highlighted the importance of having 
tsunami and geodetic data just above the source in order to better constrain the 
earthquake slip distribution 1011 . Hence, different weights are assigned to each tsunami 
and geodetic station (further details in Supplementary Information). In addition, due 
to the different behaviour of the cost functions for tsunami and geodetic data, we 
further assign different weights to both the entire data sets in order to avoid a possible 
unbalancing of the cost functions during the joint inversion. 

The average slip model is calculated as a weighted mean of a subset of the explored 
models. This subset consists of the models with the lowest cost functions (0.5% of the 
total) and which fit satisfactorily the data. The weights are inversely related to the cost 
function, thus better models count more than the others. Dispersion of the model 
parameters around their average values is assessed by performing an a-posteriori 
analysis of the explored models ensemble 10,14,16,32,71 , where the weighted standard 
deviation is taken as the error in the corresponding parameter. In this analysis the 
coefficient of variation (i.e. the standard deviation divided by the average slip value) is 
also computed for each subfault parameter (Fig. 3c). 
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